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Abstract. We survey the few exact results on the Stormer problem describing 
the dynamics of charged particles in the Earth magnetosphere. The analysis 
of this system leads to the the conclusion that charged particles are trapped 
in the Earth magnetosphere or escape to infinity, and the trapping region is 
bounded by a torus- like surface, the Van Allen inner radiation belt. In the 
trapping region, the motion of the charged particles can be periodic, quasi- 
period or chaotic. The three main effects observed in the Earth magneto- 
sphere, radiation belts, radiation aurorae and South Atlantic anomaly, are 
described in the framework described here. We discuss some new mathemat- 
ical problems suggested by the analysis of the Stormer problem. 
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1. Introduction 

Stellar and planetary magnetic environments or magnetospheres are generated 
by the motion of charged particles inside the core of stars and planets. In the 
magnetosphere of the Earth, incoming charged particles have intricate trajectories 
and are in the origin of observable radiation phenomena as is the case of radiation 
aurorae (Stormer, 1955), the Van Allen inner radiation belts (Van Allen and Frank, 
1959), and the South Atlantic anomaly (Underwood et ai, 1994). 

The magnetic field of the Earth has a strong dipolar component, (Rikitake 
and Honkura, 1985), and it is believed that, for low altitudes (< 3000 km), the 
radiation phenomena occurring in the Earth magnetosphere can be understood 
by studying the motion of nonrclativistic charged particles in a pure dipole field. 
If the dipolar component of the Earth magnetic field is considered aligned with 
the rotation axis of the Earth, the equations of motion of a charged particle in 
the dipolar field of the Earth reduce to a non-linear autonomous Hamiltonian 
dynamical system. This is the Stormer problem, (Stormer, 1955). 
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The analysis of the Stormer problem presents big challenges from the theo- 
retical, applied and computational points of view. 

From the computational point of view, the determination of the trajectories of 
high energy charged particles in the Earth magnctosphere for long periods of time 
is inaccurate and time consuming, being difficult to extract information about the 
several aspects of the radiation phenomena observed in the Earth magnetosphere. 
For example, in order to explain some of the dynamic aspects associated with 
aurorae and magnetic mirrors, adiabatic ad-hoc arguments have been introduced 
into the theory, and the motion of charged particles in the magnetosphere has been 
assumed similar to the cyclotron type motion in constant magnetic fields, (Van 
Allen and Frank, 1959), (Hess, 1968) and (Daly, 1988). 

The long lived and transient radiation belts observed in the Earth magne- 
tosphere have adverse effects on the electronics of spacecrafts, and affect com- 
munications, (Stassinopoulos and Raymond, 1988) and (Daly, 1988). Therefore, a 
qualitative and quantitative understanding of the Stormer problem has important 
applications. 

The first theoretical studies about the properties of the trajectories of charged 
particle in a dipolar field where done by DeVogelaere (1958) and Dragt (1965). In 
the work of these authors, the existence of a trapping region for charged particles 
in the dipole field of the Earth was implicitly established. 

Based on the qualitative theory of conservative maps of the plane, Dragt and 
Finn (1976) carried a comparative study between the phase space topology of the 
orbits of a generic area-preserving map of the plane and the numerically computed 
Poincare sections of the Stormer problem. They have presented numerical evidence 
about the existence of homoclinic points in the Poincare sections. According to 
these authors, this shows that the Stormer problem is insoluble, implying that 
the adiabatic magnetic moment series diverges, a basic theoretical argument used 
by Van Allen to explain some features of the radiation phenomena in planetary 
magnetospheres. The KAM approach to the Stormer problem has been developed 
by Braun in a sequence of papers, (Braun, 1970), (Braun, 1979) and (Braun, 
1981). Within this approach, it has been shown that trapped particles can have 
quasi-periodic motion, and can penetrate arbitrarily close to the dipole axis. 

Due to its intrinsic difficulty, the analysis of the Stormer problem has been 
done using a mixture of analytical and numerical techniques. Here, we are in- 
terested in surveying the exact results on the Stormer problem, separating the 
results that are numeric from the exact ones. All the exact results are summarized 
in Propositions 3.1 and 4.1. We have made extensive simulations of trajectories of 
charged particles in the Earth magnetosphere, and we have obtained the shape of 
the trapping regions for charged particles — Van Allen inner radiation belts. This 
contrast with the usual approach used in radiation environment studies, where 
radiation belt boundaries are correlated with the dipole field lines, (Daly, 1988). 

This paper is organized as follows. In the next section, we derive the equations 
of motion for the Stormer problem and we obtain its conservation laws. In section 
3, we study the motion of charged particles in the equatorial plane of the Earth. 
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In section 4, we analyse the general case for the motion on the three-dimensional 
configuration space. In the final section 5, we summarize the main results from the 
theoretical and applied points of views, and we discuss some of the mathematical 
problems suggested by the analysis done previously. 



2. Equations of motion and conservation laws 

The equation of motion of a nonrelativistic charged particle of mass m and charge 
q in a magnetic field B has the Lorentz form, 



mr 



(fx Bj . (2.1) 



We use the international system of units and B is measured in Tesla. At the surface 
of the Earth, B is in the range (0.5 x 1CT 4 - 1.0 x 10~ 4 ) Tesla = (0.5 - 1.0) Gauss. 

Magnetic fields are produced by moving charges and currents. A magnetic 
dipolc field can be produced by a current loop on a planar surface, and the resulting 
field is proportional to the current intensity times the area delimited by the current 
loop, ((Fcynman, 1964), pp. 14-7). For a current loop in the horizontal a;y-planc, 
flowing counterclockwise with current intensity /, the dipole momentum is \x z — 
fie Z} where fj, = I x area of the loop, and e z is the unit vector of the z axis. 
This current loop produces a dipole field with a dipolc momentum pointing in the 
positive direction of the z-axis. This dipolc field derives from the vector potential, 

* 11-., 11,- 1 / - -s / N 

x e r = -; x r = M z — (-ye x + xe y ) (2.2) 



4tt£oC 2 r 2 47re c 2 r 3 



where r — x 2 + y 2 + z 2 , M z is the (scalar) dipole momentum, for short, and the 
vector potential is independent of the z coordinate. As B = rotA, by (2.2), the 
Lorentz equation (2.1) describing the motion of a charge particle in a dipolc field 
is, 



3a^{yz - zy) - ay^ 



y = —3a—(xz — zx) + ax— (2-3) 
z = 3a^(xy - yx) 

where a = qM z /m. For the Earth dipolar field, the dipole momentum is M z — 
7.9 x 10 25 G cm 3 = 7.9 x 10 15 T m 3 ((Rikitake and Honkura, 1985), 1975 IGRF 
value), and, for electrons and protons, we have, 

{-1.45 x 10 27 m 3 /s (electrons) 
(2.4) 
7.88 x 10 23 m 3 /s (protons) . 

We now rescale the system of equations (2.3). With X = x/r , Y = y/r a , 
and Z = z/r a , where r = 6378136 m is the radius of the Earth, we rewrite the 
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system of equations (2.3) in the form, 



X = 3a 1 -^(YZ-ZY)-a 1 Y± 

Y = -3a 1 -^(XZ-ZX) + a 1 X^ (2.5) 
Z = 3a 1 ^(XY-YX) 



where R = ^ X 2 + Y 2 + Z 2 , 



a \ -5.588 x 10 6 s' 1 (electrons) 

a i = 4 = i , i , , ( 2 - 6 ) 



3.037 x 10 3 s" 1 (protons) . 

In this rescaled coordinate system, if R(t) = 1, the charged particle hits the surface 
of the Earth. 

The Lorentz equations (2.5) can be derived from a Lagrangian. By standard 
Lagrangian mechanics techniques, and by (2.2), we have, 



t 1 / • 2 • 2 • 2 \ ■* a* ^ i • 2 -2 • 2 \ TYlOt . . , 

i = -m(x +y + z ) + = -771(2; +y + z ) - — j- (xy - yx) 

= l - m {X 2 +Y 2 + Z 2 )-^(XY -YX). 
The conjugate momenta to the coordinates X, Y and Z are, 



(2.7) 





dL 


Y 


Px = 


dX 


— mX — max — ^ 
it J 




dL 


X 


PY = 


dY 


= mY + ma\ 




dL 


= mZ 


PZ = 




dZ 





(2.8) 



(2.9) 



and the Hamiltonian is, 

H = p.R-L = ^m(X 2 + Y 2 + Z 2 ) 

1 , 2 o 2 , X Y ma 2 X 2 + Y 2 
= ^(PX +PY +Pz)-a lPY - + a lPx - + - _. 

Hence, the system of equations (2.5) has the conservation law, 

H = ^m(X 2 + Y 2 + Z 2 ) . (2.10) 

We show now that the system of equations (2.5) has a second constant of 
motion. 

To determine the second constant of motion, we introduce cylindrical coor- 
dinates. With, X = p cos <fi, and Y = psmcfi, the Lagrangian (2.7) becomes, 

L' = \m(p 2 + p 2 ^ 2 + Z 2 ) + ^g-p 2 ^ (2.11) 
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where R = \J p 1 + Z 2 . In this coordinate system, the conjugate momenta become, 



ol 

dp 
dL 



P »=dp- = mp 



mp 2 6 + maip 2 



i? 3 



(2.12) 



dL ■ 
pz = — - = mZ 
8Z 



and the new Hamiltonian is now, 

// v-' • > ' • /:■■■■ — ' 2 • : ^ 

2 2m 



P P + Pz + \ -—-ma lR3 



As the Hamiltonian (2.13) is independent of <fr, we have p<f, = — 
(2.12), the second conservation law is, 

dL o ■ 2 1 

p^ = — r = mp 6 + maip — = c 2 = constant . 
86 Li 



dH 



. (2.13) 
0, and by 

(2.14) 



As p^ — C2 is the second equation of motion, we can write the Hamiltonian (2.13) 
in the form, 

2\ 



1 



H = - — p p + pz + mu\ 



2m 



C2 



P 



R 3 



(2.15) 



Then, the equations of motion of a nonrelativistic charged particle in a dipolc field 
reduce to, 



dH 1 



dH 1 



^ dp p m Pp 



Pp 
Z = 
Pz 



2 2 p° 2 P P 

c 2 + Sma^ - ma x - 3aic 2 -^ 



dp mp 3 
dH 1 

a — = —P z 

dpz rn 

= -|| =3ma?^#-3aic 2 z 



or, 



1 



P - 



2 , o 2 P 2 P 



, «lC2 P 

m R 5 



3a 



p Z Ct\C2 Z 



(2.16) 



1 R s m R 5 
together with the conservation laws, 

1 / o o 



2m? 



C2 P 
Pp" +PZ~ + \ — - ma 1^3 

1 



Cl 



(2.17) 



mp 2 ^ + moi\p 2 -^ = c 2 
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where R = \J p 1 + Z" 2 . Note that, we have rescaled the Hamiltonian of the Stormer 
problem to, 

H eff = \ + 2?) +^(j-m*^) 2 =T + V eff ( P , Z) (2.18) 

where T is a scaled kinetic energy. 

Hence, the motion of charged particle in a dipole field is described by a two- 
degree of freedom Hamiltonian system. The effective Hamiltonian is parameterized 
by the second conservation law in (2.17), which determines the time dependence 
of the angular cylindrical coordinate. 

Before analysing the general topology of the orbits of the phase space flow 
of equations (2.16), we first consider the case where the motion restricted to the 
plane Z = 0. 



3. Motion in the equatorial plane of the dipole field 

Here, we consider the case where charged particles are constrained to the equatorial 
plane of the Earth, the plane Z = 0, for every t > 0. In this case, by (2.16), the 
equations of motion reduce to, 



£2|i_ _ 3 ^i£20 1 + 2a2 1 = _dVe £l 
m z p A m p 4 p° dp 



where C20 is the value of the constant C2 evaluated at Z = 0, and, by (2.18), the 
potential function V e //(p) is, 

Introducing pz = and Z = into (2.17), the two conservation laws reduce to, 

2\ 



Pp 2 + I — - ma i^ I I = C10 



C20 



2m2 I" ' V P V 



1 

mp <p + m«i - = c 2 o 



(3.3) 



P 



where, 



1 



Cw = ^ 0)2 + ^{pfj- mai 7(oj 

2 • 1 

c 2 o = mp(0) 0(0) + mar 



( _£H2_ 



(3.4) 



P(0) 

and C10 is the value of the effective total energy evaluated at Z = 0. 

Integrating the second conservation law in (3.3) by quadratures, we obtain, 
for the angular coordinate 4>(t), 

^= m+ lX^-w?) dt - (3 - 5) 
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Figure 1. Effective potential associated with the motion of a 
charged particle (protons, a\ > 0) in the equatorial plane of the 
Earth (Z = 0). a) If c 2 q > and p > 0, the effective potential has 
a maximum for p = p 2 , a minimum for p — pi, V e ff(pi) — 0, and 
Veffipi) = 4 /(32™ 4 a 2)- At p = Pa = (V2 - 1) P2 , V eff (p ) = 
V e ff(p 2 ). b) If C20 < and p > 0, V e ff(p) is a monotonically 
decreasing function of the argument, and the effective energy sur- 
faces have no compact components. 



So, if the solution p(t) of equation (3.1) is known, the temporal dependency of 
angular coordinate is obtained from (3.5). 

In these conditions, the motion of the charged particle in the equatorial plane 
of the Earth is completely determined by equation (3.1), derived from the effective 
Hamiltonian, 

tt / • \ rp.fr 1 -2 , C 2Q 1 QlC 2 o 1 ttjl 

Heff(p,p)=T + V eff = -p +—---—- + -- (3.6) 

where T is a scaled kinetic energy. If C20 > and for p > 0, the potential function 
V e ff(p) has one local minimum and one local maximum at, 

Pi = and p2 = 2 (3.7) 

C20 c 20 

respectively, Figure la). By direct calculation, we have, V e ff(p2) = c|o/(32™ 4 a^), 
and Veff(pi) = 0. Therefore, for < H e ff(p,p) < V e ff(p2), the constant effective 
energy surface contains a compact component. 

If C20 < and for positive values of p, the potential function V e ff (p) is a 
monotonically decreasing function of the argument, Figure lb). Then, we have: 

Proposition 3.1. In the equatorial plane of a dipole field, a nonrelativistic (posi- 
tively) charged particle precesses around the dipole axis, along a circle with radius 
pi, provided: (i) H eff (p(0),p(0)) < V e ff(p 2 ), c 20 > 0, p(0) ^ p\ and p(0) < p 2 , 
where H e ff is defined in (2>.Q), the constants p\ and p 2 are given in (3.7), and the 
constant c 2 o is defined in (3A). For initial conditions near the circumference of 
radius p = p\, the Larmor or precession period is, Tj, = 27rm 3 af /c 20 - The phase 
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advance per Larmor period is, 

and the period of rotation around the Earth is T r = 2ttTl / A<j>. 

(ii) If H eff (p(0),p(0)) = V eff (p 2 ), c 20 > 0, andp(0) < p 2 , or, p(0) > p 2 and 
p(0) < 0, then, p(t) — > p 2 , as t — > 00. 

(Hi) If p(0) = p 2 [resp., p(0) = p\], p(0) = 0, (f)(0) ^ 0, and c 2 o > 0, then the 
charged particle has a circular trajectory with radius p = p 2 [resp., p = p\], and 
the period of rotation around the Earth is T r — 2ir/<p(0). If (f)(0) = 0, the charged 
particle is at rest. 

(iv) If H ej f(p(0),p(0)) > V ef f(p 2 ) andc 20 > 0, or, c 20 < 0, then, p(t) — > 00, 
as t — > 00, and we have escape trajectories. 

Proof. To prove the proposition, we must ensure first that the solution of the 
differential equation (3.1) exists, and is defined for every t > 0. In the cases (i)- 
(iii), the existence of solutions for every t > follows because the initial conditions 
on phase space are on the compact components of the level sets of the Hamiltonian 
function (3.6), (Chillingworth (1976), pp. 187; Nemitskii and Stepanov (1960), pp. 
8). 

A simple phase space analysis shows that the vector field associated with 
equation (3.1) has a centre type fixed point with phase space (cylindrical) coor- 
dinates (pi,p = 0). This centre type fixed point is inside the homoclinic loop of 
a saddle point with coordinates (p 2 ,p = 0). The conditions in (i) correspond to 
initial conditions inside the homoclinic loop, and away from the centre fixed point, 
p(0) 7^ p\. To calculate the Larmor frequency and the phase advance per Larmor 
period, we linearize equation (3.1) around pi, and we obtain, 

c 6 

x+^ 1 x = 



where x = p— pi, and 0%/ '(m 6 af) = lo\. The Larmor period is Tl = 2-k/lol- Note 

d 2 V e ff 

that, by a direct calculation, u\ = - & 2 — (p = pi). 

By (3.5), the phase advance per Larmor period is, 



A< ^ So (mp(s) 2 p(s) 3 ) dt 

From the above linearized differential equation, in the vicinity of the circumfer- 
ence of radius p = pi, p(t) = pi + (p — pi) cos(wLt) + p(0) sin(w ^t) / 'uj l and after 
substitution into the expression for A<j), we obtain the result. The rotation period 
around the Earth is obtained from the conditions, nA(f> = 2ir and T r — nT^. 

The conditions in (ii) correspond to initial conditions on the stable branches 
of the homoclinic orbit of the saddle point with coordinates (p 2 ,p~ 0). 
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For case (Hi) , as (p2,p = 0) is a fixed point of the differential equation (3.1), 
if <^(0) 7^ 0, the charged particle has a circular trajectory around the origin of 
coordinates. To calculate the period of the trajectory, by (3.5), we have, 

\fnpi p%) 

Introducing the value of the constant c 2 o into the above expression, by (3.4), we 
obtain T r — 2n/<j>(0). For p(0) = pi, the proof is similar. 

In case (iv), the effective potential function monotonically decreases as p 
increases from p = 0, and we are in the conditions of escape trajectories. The pro- 
longation of solutions for every t > 0, follows from the condition that fi(p)/\\p\ | — > 
constant as p — > oo (Ncmitskii and Stepanov (1960), pp. 9). □ 

By Proposition 3.1, the trajectory of a charged particle on the equatorial 
plane of an axially symmetric magnetosphere can be unbounded, can be trapped 
in an annular region around the Earth, or can collide with the surface of the Earth. 
As po, pi and p 2 depend on c 2 o, there are charged particles with different energies 
trapped in the annular region [po,p 2 \. This annular region is the equatorial cross 
section of the Van Allen inner radiation belt, (Van Allen and Frank, 1959). 

As po, pi and p 2 depend on the initial conditions p(0) and </>(0) through c 2 o, 
for the same belt parameter, there are particles that escape from the trapping 
region, escaping to infinity or hitting the surface of the Earth. This is due to the 
fact that the belt parameters are independent of p(0). By Proposition 3.1- (iv), and 
by the conservation of energy, a charged particle escapes from the inner Van Allen 
radiation belt if, 

|P(0)| > ^2\V eff ( P2 )~V eff (p(0))\ (3.8) 

By the same argument, if |p(0)| > \c 2n /m — ai\, a particle with initial velocity 
p(0) < at infinity, and velocity vector within the equatorial plane of the Earth, 
after an infinite time, hits the surface of the Earth. 

In Figure 2, we show the limits of the equatorial cross section of a Van Allen 
inner radiation belt, the trajectories of trapped protons, and an escape trajectory. 
In the three cases shown, the Van Allen parameters are the same. The trajectories 
shown in Figure 2 have been calculated by the numerical integration of equation 
(3.1), (see the Appendix). The angular variable <f) has been obtained from the 
discretization of (3.5), 

<f> n+1 = 4> n + At - ^f) (3.9) 

where At is the discretization time step, <p n = (p(nAt), p n = p(nAt), and n — 
0,1,- ■■■ 

As the radiation belt parameters pa, pi and p 2 depend on the angular velocity 
4>(0), we can have several Van Allen radiation belts at different altitudes, and 
trapped particles with different energies. 
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Figure 2. Equatorial cross section of the Van Allen inner ra- 
diation belt, for protons in the equatorial plane of the Earth. 
We show the trajectories of three protons with initial conditions 
(cylindrical coordinates): a) p(0) = 3.0, p(0) = 10.0, 0(0) = 0.0 
and 0(0) = 10.0; b) p(0) = 3.0, p(0) = 80.0, 0(0) = 3tt/2 
and 0(0) = 10.0; c) p(0) = 3.0, p(0) = 100.0, 0(0) = tt and 
0(0) = 10.0. The first and the second trajectories, a) and b), 
correspond to the precession of protons around the Earth. The 
third case is an escape trajectory. By (3.8), the escape condition 
is |p(0)| > 95.42. For the three trajectories shown, the radiation 
belt parameters are po — 2.282, p\ = 2.755, p 2 = 5.510, and 
c 2 o = 1.844 x 10~ 24 . The total effective energies are: a) Cio = 500; 
b) Cio = 3650, and, c) c w — 5450. The trajectories have been 
calculate with the Stormer-Verlet numerical method (Appendix) 
with the time step At = 0.0001, and the angular coordinate has 
been calculated by (3.9). For protons with precessing trajectories 
near the circumference of radius pi, by Proposition 3.1, the Lar- 
mor period is Tl = 0.043 s. The first two particle trajectories have 
been calculated from t = up to t = 2 s. 
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4. Motion in the three-dimensional space 

The equations of motion of a charged particle in a dipole field — equations (2.16), 
are derived from the effective potential function, 

I m.rv^ 



where R = \J p 2 + Z 2 . In Figure 3a), we show the graph of the potential function 
V e ff(p,Z), and its level lines, for c 2 > 0. 

As we have seen in the previous section, on the equatorial plane Z = 0, the 
exterior boundary of the trapping region is the circumference of radius p 2 , which 
corresponds to the unique local maximum of the potential function V e ff(p). The 
equation V e ff (p, Z) = V e ff(p 2 ) defines a bounded region in the (p, Z) plane whose 
closure is a compact set of maximal area 1 , Figure 3b). The compact components 
of the level sets of the effective Hamiltonian function (2.18) are obtained as the 
topological product of a compact plane set with the compact sets in the interior of 
the region defined by the equation V e f f(p, Z) = V e ff(p 2 ). Therefore, any charged 
particle with an effective energy on a compact component of the effective Hamil- 
tonian is trapped in a torus-like region around the Earth. This trapping region is 
the Van Allen inner radiation belt of the Earth. 

In the rescaled coordinate system introduced in Section 2, the level lines 
V e ff(p, Z) — V e ff(p 2 ) hit the surface of the Earth for yj p 2 + Z 2 = 1. So, we 
define the contact points of a specific radiation belt as the points at the surface 
of the Earth where V e ff(p,Z) = 0. If c 2 > 0, this potential function takes its 

c 2 P 

minimum value (V e fj = 0) along the line — = mai — , Figure 3b), and, for 
R = 1, we have, 



Pc = x \— = ^= (4.2) 
mai ^/Pi 



which corresponds to the latitudes, 



9 C = arccospc = arccos . / ° 2 . (4.3) 

We now summarize the main features of the dynamics associated with the 
two degrees of freedom Hamiltonian (2.18). 

Proposition 4.1. We consider the motion of a nonrelativistic (positively) charged 
particle in a dipole field. If one of the coordinates Z(0) or Z(0) of the initial condi- 
tion is different from zero, and (p(0), Z(0)) ^ (0,0), then the charged particle has 
bounded motion, being trapped in a torus-like region around the Earth, provided: 



lr The compact set is obtained by adding the point (0, 0) to the open set defined by V e ff(p, Z) = 
V e ff(p2)- In the following, compact level sets are always obtained by adding this exceptional 
point to the open and bounded level sets of the effective potential function V e //(p, Z). 
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Figure 3. a) Graph of the effective potential function (4.1), with 
C2 > 0, for the three dimensional Stormer problem for protons, 
b) Effective potential energy level line (heavy lines), V e ff(p, Z) = 
Veff(P2)- For particles characterized by the constant of motion 
c 2 > 0, the Van Allen inner radiation belt is the toric-like surface 
obtained by rotating the level lines shown in b) around the axis of 
the dipole field of the Earth. The annular region shown in Figure 
2, delimited by the circumferences of radius po and p 2l is the 
Z = cross section of the Van Allen inner radiation belt. The 
dotted line is the local minimum of the effective potential energy 
function V e f / (p, Z) . 

(i) H eff (p(0),p(0),Z(0),Z(0)) < V eff (p 2 ), c 2 > 0, andp(O) < p 2 , where p 2 is de- 
fined in (3.7), the Hamiltonian function H e ff is defined in (2.18), and the constant 
of motion c 2 is defined in (2.17). In the particular case where, Z(0) = 0, p(0) = 0, 

and c 2 /p(0) — ma\p(fS) / ' \J p(0) 2 + Z(Q) 2 , the charged particle remains, for every 
t > 0, m the plane Z = Z(0), and Z(0) G [-2ma 1 /(3V3c 2 ), 2mai/(3v / 3c 2 )] . If 
<j)(0) 7^ 0, it rotates around the dipole axis. If 0(0) = 7 the particle is at rest in 
the plane Z = Z(0). 

(ii) If one of the coordinates Z(0) or Z(0) of the initial condition is different from 
zero, and one of the above conditions is not verified, then p(t) — > oo, as t — > oo, 
and we have an escape trajectory. 

Proof. If Z(0) = Z(Q) = 0, the motion is restricted to the plane perpendicu- 
lar to the dipole axis, and we obtain the case of Proposition 3.1 of the previous 
section. If (p(0), Z(0)) = (0,0), the equation of motion (2.16) is not defined. In 
case (i), the initial condition is in a compact component of the effective Hamil- 
tonian, and the motion is bounded for every t > 0. If, Z(0) = 0, p(0) = 0, 

and c 2 /p(0) = maip(0)/\/ p(0) 2 + Z(0) 2 , where c 2 is defined in (2.17), and 
Z(0) £ [— 2moei/(3\/3c 2 ),2mai/(3\/3c 2 )], the differential equation (2.16) has one 
or two fixed points in the plane Z = Z(0). 
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In case (ii), the level sets of the Hamiltonian function are not bounded, and 
the prolongation of solutions as t — > oo follows as in case (iv) of Proposition 



In Figure 4, we show the trajectories of several particles trapped in a Van 
Allen inner radiation belt characterized by the constant ci > 0. The trajectories 
have been calculate with the Stormer- Verlet numerical method (Appendix) with 
integration time step At — 0.00002, and total integration time t = 2s. From (2.17), 
it follows that the angular coordinate <p(t) is given by, 



where At is the discretization time step, <f) n = 4>{nAt), p n = p(nAt), Z n = Z(nAt), 
and n = 0, 1, • • • . In the trajectories of Figure 4, the angular coordinate 4>{t) has 
been calculated with (4.4). 

In the plots on the right-hand side of Figure 4, we show the projection of 
the boundary of the level set of the effective Hamiltonian on the (p, Z) plane. 
The trajectories of trapped charged particles are in the interior of these boundary 
curves. Particles with low effective energy have trajectories concentrated in the 
vicinity of the equatorial plane of the Earth, Figure 4a). Increasing the effective 
energy of the particles, their trajectories approach the polar regions of the Earth 
and eventually hit the surface of the Earth. By (4.3), and for the parameter values 
of Figure 4, the contact points of the Van Allen inner radiation belt with the 
Earth are located at the latitudes 9 C = ±53.76°. These simulations suggest that 
the constant effective energy surface is not filled density by a unique trajectory. 

To characterize more precisely the dynamics of the charged particles trapped 
in the dipole field of the Earth, we can eventually use KAM techniques, (Licht- 
enberg and Lieberman, 1983), or construct a Poincare map for the equations of 
motion (2.16). 

To pursue a KAM approach, (Lichtenberg and Lieberman, 1983), the first 
step is to find an integrable Hamiltonian system leaving invariant a two-dimensional 
torus in the four-dimensional phase space of the differential equations (2.16). 
For that, we develop in Taylor series around the points (p = pi,Z = 0) and 
(p = p2,Z = 0) the second members of the differential equations in (2.16). By 
Proposition 3.1, at these points the motion is integrable and periodic in the equa- 
torial plane of the Earth. In the first case, the differential equations (2.16) become, 



3.1. 



□ 




and by discretization, we obtain, 




(4.4) 




(4.5) 
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Figure 4. Trajectories of three protons in the Van Allen inner 
radiation belt, with initial conditions: a) p(0) = 3.0, p(0) = 10.0, 
0(0) = 0.0, 0(0) = 10.0, Z(0) = 0.5 and Z(0) = 0.0; b) p(0) = 3.0, 
p(0) = 30.0, 0(0) = 0.0, 0(0) = 10.0, Z(0) = 0.5 and Z(0) = 0.0; 
c) p(0) = 3.0, p(0) = 80.0, 0(0) = 0.0, 0(0) = 10.0, Z(0) = 0.5 
and Z(0) = 0.0. The radiation belt parameters are po — 2.370, 
pi = 2.861, p 2 = 5.722, and c 2 = 1.776 x 10~ 24 The total effective 
energies are: a) c\ = 500; b) c\ = 900, and, c) ci = 3650. 
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and, in the second case, we obtain, 



P 



32m 6 af 
3c§ 



(P- Pi) 



< 



(4.6) 



Z 



64m 6 af 



Z. 



As the solutions of the linear equations (4.5) and (4.6) are unbounded, we loose the 
property of boundness already contained in Proposition 4.1- (i) and also the possi- 
bility of having a family of invariant two-dimensional torus in the four-dimensional 
phase space of the unperturbed systems (4.5) and (4.6). 

The other way of characterizing the dynamics of the charged particles in the 
Van Allen inner belt is to find a Poincare map for the equations of motion (2.16). 

We consider that at most one of the coordinates Z(0) and Z(0) of the initial 
condition of the differential equation (2.16) is different from zero, and the particle 
has bounded motion (c 2 > 0). To avoid degenerate situations, we also consider 
that, (p(0), Z(Q)) 7^ (0,0) and the effective energy function is positive, c\ > 0. 
Under these conditions, by Proposition i.l-(i), the initial condition is on a compact 
component of the level sets of the Hamiltonian function (2.18). So, the solution of 
the differential equation (2.16) exists, and is defined for all t > (Chillingworth 
(1976), pp. 187; Nemitskii and Stepanov (1960), pp. 8). 

The compact components of the level sets of the Hamiltonian function H e f / 
are obtained as the topological product of a compact plane set with the com- 
pact level set of the potential function V e ff(p, Z). These level sets are compact, 
provided C2 > 0. The compact level sets of the Hamiltonian function are three- 
dimensional compact manifolds embedded in the four dimensional phase space, 
and their intersection with the three-dimensional hypcrplanc Z = is a two- 
dimensional compact manifold X C2 . This two-dimensional manifold is a compact 
set in the four dimensional phase space of the differential equation (2.16) and, by 
(2.18), has local coordinates p and p. 

In the conditions of Proposition A.l-(i), any orbit initiated in an initial condi- 
tion on the effective energy level sets can cross the hyperplane Z = 0, intersecting 
transversally the two-dimensional compact manifold X C2 . Therefore, the compact 
two-dimensional manifold E C2 is a good candidate for a Poincare section of the 
differential equation (2.16). This construction is used to justify the existence of 
two-dimensional conservative Poincare maps for two-degrees of freedom Hamil- 
tonian systems, (Lichtenberg and Lieberman (1983), pp. 17-20). However, there 
arc two main difficulties in proving that £ C2 is the domain of a Poincare map. 
The first case, it is difficult to prove that any trajectory that crosses transversally 
the section S C2 will return to it after a finite time. The second difficult questions 
concerns the unicity of the trajectories crossing the plane Z = 0. 

From Proposition 4.1- (i), it follows that there are orbits of trapped particles 
that never cross E C2 . However, due to the invariance of the differential equation 
(2.16) for transformations of t into —t, any initial condition on the plane Z = 
with Z 7^ 0, has a prolongation for positive and negative values of t, crossing 
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Figure 5. Poincare maps for the constant of motion C2 = 1.776 x 
1(T 24 . The Poincare map have been calculated with the initial 
conditions: a) p(0) = 3.0, 0(0) = 0.0, 0(0) = 5.47089, Z(0) = 0.0 
and Z(0) = 0.5; b) p(0) = 3.0, 0(0) = 0.0, 0(0) = 10.0, Z(0) = 0.5 
and Z(Q) = 0.5. In both figures, the initial conditions p(0) have 
been changed in the interval [0,85]. The dots have coordinates 
(pi,0) and (p2)0), with pi < p 2 - Both Poincare sections have 
been generated with different families of initial conditions. This 
suggests that there are different particle trajectories with the same 
effective energy that cross the two-dimensional section £ C2 at the 
same point. 



transversally the plane Z = at most once. With this property, we can test 
numerically the structure of the orbits that cross the set £ C2 . In Figure 5, we 
show the computed orbits on the Poincare section £ C2 , for two families of initial 
conditions and the same value of the constant of motion c 2 . 

In Figure 5a), we have chosen initial conditions on the plane Z(0) = 0, with 
Z(0) 0. For all the analysed initial conditions, the orbits crossed several times 
the surface of section E C2 . Due to the structure of orbits, the motion appears to be 
quasi-periodic. However, if the initial conditions are away from the plane Z(0) = 
and for the same value of C2, the topology of the orbits on the section S C2 changes 
drastically, Figure 5b). Comparing Figure 5a) with 5b), the structure of the orbits 
arc incompatible, in the sense that different orbits cross S C2 at the same point. 
This shows that S C2 can not be the domain of a Poincare map for the differential 
equations (2.16). On the other hand, the structure of orbits in Figure 5b) suggests 
the existence of transversal homoclinic intersections and therefore chaotic motion 
inside the trapping region of the Stormer problem. 
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5. Conclusions 

The Hamiltonian dynamical system describing the motion of a charged particle 
in a dipole field (Stormer problem) has been reduced to a two-degrees of freedom 
system. This reduced system has two constants of motion. It has been shown that 
the trajectories of charged particles can be periodic and quasi-periodic, and that, 
for a suitable choice of the initial conditions in between a minimal and maximal 
height from the equatorial plane, the planes perpendicular to the dipole axis are 
invariant for the motion of charged particles. 

From a more global point of view, the trajectories of the charged particles in 
the Earth dipole field can be trapped in a torus-like region surrounding the Earth, 
or can be scattered and escape to infinity. The torus-like trapping region around 
the Earth can be interpreted as the Van Allen inner radiation belt, as measured 
by particle detectors in spacecrafts. The physical effects associated with radiation 
phenomena intrinsic to accelerated charged particles (Bremsstrahlung) account for 
the phenomena of radiation aurorac. The numerically computed trajectories of 
trapped charged particles suggests the existence of chaotic motion inside the Van 
Allen inner belts around the dipole axis of the Earth. The properties of the Stormer 
dynamical systems are summarized in Propositions 3.1 and 4.1. 

From the applied physics point of view, the Stormer problem deviates from 
the real situation in essentially three ways. In the first case, the magnetic field 
of the Earth has a strong quadrupolar component that was not considered. The 
second drastic simplification has to do with the fact that the dipole axis is not 
coincident with the rotation axis of the Earth. The rotation of the Earth introduces 
a time periodic forcing into the equations of motion. In the third simplification, we 
have considered implicitly that charged particles do not radiate (Bremsstrahlung) 
when subject to accelerating forces. However, the three main effects observed in 
the Earth magnetosphere, radiation belts, radiation aurorae and South Atlantic 
anomaly, are described by the simplified model. 

The existence of periodic, quasi-periodic and chaotic trajectories in the Stor- 
mer problem together with the Bremsstrahlung effect, shows that radiation belts 
make a shield protection for the high energy charged particles arriving at Earth. 
These particles loose kinetic energy by Bremsstrahlung when they are scattered 
and when they are trapped in the Van Allen inner radiation belts. 

From the mathematical point of view, the Stormer problem poses some open 
problems for the dynamics of two-degrees of freedom Hamiltonian systems. The 
Stormer dynamical systems is not generated from the perturbation of a two- 
dimensional torus in a four-dimensional phase space, suggesting the existence of 
non KAM mechanism for the generation of bounded motion in two-degrees of 
freedom Hamiltonian systems. On the other hand, the arguments used in the con- 
struction of Poincare sections in two-degrees of freedom Hamiltonian systems with 
two conservation laws lead to the non-unicity of the orbits on these surfaces of 
section. 
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Appendix 

To integrate numerically the equations of motion (2.16) and (3.1), we have used 
the explicit Stormer-Verlet method of order 2, (Haircr et al. (2002), pp. 14 and 
177), 

Pn+l/2 =Pn~ M—{qi) 

a 9H. 

Qn+1 =Qn+ At—{ Pn+1/2 ) 

A^fe +1 ) 



Pn+l — Pn+l/2 



where H(qi, . . . , q n ,Pi, ■ ■ ■ ,Pn) is the Hamiltonian function, and At is the integra- 
tion time step. For example, for equation (3.1), the Stormer-Verlet method reduces 
to, 

At 2 

Pn+l = Pn + Atp ri + j '{fin) 

Pn+l = Pn + -y/lPn) + "y/(Pn+l) 
dV e ff 

where p n = p(nAt), p n = p(nAt), f(p) = ^LL. xhis method is symplectic, 

dp 

being area preserving. 

The Stormer-Verlet method has the advantage of being explicit, and the 
integration accuracy is obtained by decreasing the time step At. However, it does 
not conserve the energy function. We have tested other higher order numerical 
integrators but, with this Stormer-Verlet method, the overall behaviour of the 
solutions is closer to the exact results. For a detailed theoretical discussion about 
the Stormer-Verlet method see Tupper (2005). 
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